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Abstract 

The objective of this effort is to develop an efficient and accurate computational 
methodology to predict both detailed thermo-fluid environments and global 
characteristics of the internal ballistics for a hypothetical solid-core nuclear thermal 
thrust chamber assembly (NTTCA). Several numerical and multi-physics thermo- 
fluid models, such as real fluid, chemically reacting, turbulence, conjugate heat 
transfer, porosity, and power generation, were incorporated into an unstructured- 
grid, pressure-based computational fluid dynamics solver as the underlying 
computational methodology. The numerical simulations of detailed thermo-fluid 
environment of a single flow element provide a mechanism to estimate the 
thermal stress and possible occurrence of the mid-section corrosion of the solid 
core. In addition, the numerical results of the detailed simulation were employed 
to fine tune the porosity model mimic the pressure drop and thermal load of the 
coolant flow through a single flow element. The use of the tuned porosity model 
enables an efficient simulation of the entire NTTCA system, and evaluating its 
performance during the design cycle. 

I. Introduction 

Nuclear thermal propulsion can carry far larger payloads and reduce travel time for 
astronauts traveling to Mars and other planets of the solar system than any chemical propulsion 
currently available. One of the concepts that was extensively tested during the Rover/NERNA 
era and appear to be the most feasible is the solid-core concept . 1 This concept involves a solid- 
core reactor consisting of hundreds of heat generating solid flow elements, and each flow 
element containing tens of flow channels through which the working fluid acquires energy and 
expands in a high expansion nozzle to generate thrust. Hydrogen is most often chosen as the 
propellant due to its low molecular weight. The solid-core reactor therefore resembles a heat 
exchanger. To minimize the effect of its weight, the reactor often operates at very high 
temperature and power density, which imposes real challenges to the integrity of the flow 
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element material. To make the solid-core reactor a viable concept for Mars missions, we must 
understand the effect of hydrogen as a high temperature working fluid and develop materials that 
withstand the harsh flow element environment. 

One of the impacts of operating at the combination of high temperature and high power 
density is a phenomenon known as the mid-section corrosion, as reported during the legacy 
engine tests. 1, 2 There was an excessive mass loss of the flow element material near the mid- 
section during testing. The symptom was cracked coating layer while the purpose of coating was 
to isolate the carbonaceous compound in the flow element matrix from the attack by hydrogen. 
The causes of mid-section corrosion were speculated as a mismatch in the thermal expansion of 
flow element and its coating material, large temperature gradients in the solid core, and change 
of solid thermal property due to irradiation. 4 Note that another speculation was that the flow was 
choked in the long flow channels of the flow element. 

One way to help understanding the effect of hydrogen as a working fluid and developing 
materials that withstand the harsh environment is to develop computational methodology that 
can accurately predict thermo-fluid environments inside the nuclear thermal thrust chamber 
assembly and reproduce the flow element thermal environment occurring in the legacy engine 
tests. The objective of this effort is therefore to develop an efficient and accurate multiphysics 
thermal-fluid computational methodology to predict environments for a hypothetical solid-core 
thrust chamber and the associated flow element, similar to those in the Small Engine. 1 The Small 
Engine was a paper engine designed near the end of the Rover/NERVA era and bears common 
features of other legacy engines tested during that time period, but was never built nor tested. 
The hypothetical thrust chamber and flow element, which are being redesigned by the System 
Analysis Group at Marshall Space Flight Center, consists of 564 
flow elements and of 19 flow channels for each flow element, 
respectively. The computational methodology was based on an 
existing Unstructured-grid Navier- Stokes Internal-external 

computational fluid dynamics Code (UNIC 5 ' 7 ). Conjugate heat 
transfer formulations for coupling fluid dynamics and conductive 
heat transfer in solids and for flow and conductive heat transfer in 
porous media were developed and tested. The UNIC code has been 
well validated and employed to simulate a great variety of 
engineering problems ranging from internal to external flows, 
incompressible to compressible flows, single-phase to multi-phase 
flow, and inert to reacting flows. 

A two-pronged approach was employed in this effort: a 
detailed analysis of a NERVA-type 19-channel flow element, and a 
global analysis of the entire thrust chamber. Two groups of detailed 
analyses were conducted for the 19-channel flow element, as shown 
in Figure 1. The first group simulated the full-length 19-channel 
flow element with three different power generation distributions to 
investigate the occurrence of mid-section corrosion problem and 
potential flow chocking in the flow channel. These three 
distributions include pure sine function in the axial direction with Figurel: Geometry of a 
pure cosine function in the radial direction, pure sine function in the 1 9-channel flow element 
axial direction with clipped cosine function in the radial direction, 

and clipped sine function in the axial direction with clipped cosine function in the radial 



direction. The second group simulated a one-eighth-length 19-channel flow element with 
different side-wall temperature and flow channel diameters to obtain data for calibrating the 
porosity model used in the global analysis of the entire thrust chamber. In the global thrust 
chamber analysis, the 19-channels for each flow element were lumped together as a porous 
media to save computational resources, and the core surrounding components such as the slats 
and reflector were treated as heat conducting solids to provide accurate boundary condition for 
the solid-core environment. The porosity model is employed to represent the effect of the 
friction loss and heat transfer as the working fluid flowing through the flow element. The use of 
the porosity model enables an efficient simulation of the entire thrust chamber, and evaluation of 
its performance during the design cycle. 

In order to support both the detailed and global analyses, several physical submodels were 
either improved or added into the UNIC code 8-9 . These code improvements include: an anisotropic 
drag and heat transfer porosity model to account for the directional effect of the flow channel, a 
source term in the energy to model the power generation of the solid core with a user’s specified 
distributions, validation of the conjugate heat transfer capability built in the UNIC code with a well 
known thermal analysis code, SINDA 10 to provide an accurate calculation of heat transfer through 
the solid-fluid interface. The tuned porosity model for the flow element and the numerical results of 
both the detailed and global analyses are reported herein. . 


II. Numerical Methodology 
A. Computational Fluid Dynamics 

The employed CFD solver, UNIC, solves a set of Reynolds-averaged governing 
equations (continuity, Navier-Stokes, energy, species mass fraction, etc.) to satisfy the 
conservation laws for a turbulent flow of interest. The set of governing equation can be written 
in Cartesian tensor form: 
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where p is the fluid density p is the pressure, Vj - ( u , v, w) stands for the velocity components in 
x-, y-, and z-coordinates respectively, h, and h are the total and static enthalpies, k is the 
turbulence kinetic energy, n and s are the production and dissipation rates of turbulence, a, and 
Si are the mass fraction and production/destruction rate of /-th species, Q r is the radiative heat 
flux, Sv and 5/, are the source/sink terms of the momentum and energy equations, p and p, are the 
fluid and eddy viscosity, zp and z, are the viscous stress and Reynolds stress, Pr and Pr, are the 
Prandtl and turbulent Prandtl numbers, Sc and Sc, are Schmidt and turbulent Schmidt numbers, 
Ci, C 2 , C 3 , cTk, and <7 S are turbulence modeling constants. Detailed expressions for the k-e 
models and wall functions can be found in [1 1]. An extended k-e turbulence model 12 was used to 
close the system of Reynolds-averaged governing equations. A modified wall function 
approach 13, 14 was employed to provide wall boundary layer solutions that are less sensitive to 
the near-wall grid spacing. 

A predictor and multi-corrector pressure-based solution algorithm 15, 16 was employed in 
the UNIC code to couple the set of governing equations such that both compressible and 
incompressible flows can be solved in a unified framework without using ad-hoc artificial 
compressibility and/or a pre-conditioning method. The employed predictor-corrector solution 
method 5 is based on modified pressure-velocity coupling approach of the SIMPLE-type 16 
algorithm which includes the compressibility effects and is applicable to flows at all speeds. In 
order to handle problems with complex geometries, the UNIC code employs a cell-centered 
unstructured finite volume method 6, 7 to solve for the governing equations in the curvilinear 
coordinates, in which the primary variables are the Cartesian velocity components, pressure, 
total enthalpy, turbulence kinetic energy, turbulence dissipation and mass fractions of chemical 
species. 

The inviscid flux is evaluated through the values at the upwind cell and a linear 
reconstruction procedure to achieve second order accuracy. A multi-dimensional linear 
reconstruction approach by Barth and Jespersen 17 was used in the cell reconstruction to achieve 
higher-order accuracy for the convection terms. A second-order central-difference scheme was 
employed to discretize the diffusion fluxes and source terms. A dual-time sub-iteration method 
is employed for time-accurate time-marching computations. A pressure damping term, Rhie and 
Chow 18 , is applied to the evaluation of mass flux at the cell interface to avoid the even-odd 
decoupling of velocity and pressure fields. All the discretized governing equations are solved 
using the preconditioned Bi-CGSTAB 19 matrix solver, except the pressure-correction equation 
which has an option to be solved using GMRES 20 matrix solver when the matrix is ill- 
conditioned. An algebraic multi-grid (AMG) solver 21 was included such that users can activate 
to improve the convergence if desired. The UNIC code also has mesh refinement and coarsening 
capabilities to improve the numerical accuracy and computational efficiency. In order to 
efficiently simulate problems involving large numbers of meshes, the UNIC code employed 
parallel computing with domain decomposition, where exchange of data between processors is 
done by using MPI 22 . Domain decomposition (partitioning the computational domain into 
several sub-domains handled by different computer processors) can be accomplished by using 
METIS 23 or a native partitioning routine in the UNIC code. 

In addition, numerous submodels were incorporated in the UNIC code to model various 
physics and transport phenomena. These submodels include: 1) finite-rate and equilibrium 
chemistry for combustion, 2 ) conjugate heat transfer for coupling heat convection and 
conduction, 3) both the heterogeneous spray model with Eulerian/Lagrangian approach and the 
homogeneous spray model with the Eulerian/Eulerian approach and real-fluid property models 



for liquid spray flow, 4) the finite volume method for solving radiation transport equation (RTE), 
5) solving translational and vibrational energy equations for thermal non-equilibrium effect, 6) 
porosity model for flow through a porous media, etc. During the last decade, the UNIC code has 
been improved to equip models appropriate to resolve the complex physics, and has been 
employed to simulate a great variety of engineering problems. The results have shown that the 
UNIC code is as computationally efficient as any CFD codes. Hence, we employed the UNIC 
code and modified it to account for some specific effects of interest in the present study, such as 
power generation, non-isotropic porosity with or without heat conduction, variable thermal 
conductivities and specific heat based on different solid materials and temperatures. Some of 
these code improvements performed in this study are described in the following. 


B. Conjugate Heat Transfer: 

The framework of conjugate heat transfer is to solve the heat transfer in the fluid flow 
and the heat conduction in the solid in a coupled manner. The governing equation describing the 
heat transfer in the fluid flow is shown in Eq. (3), where the heat conduction in the solid can be 
written as: 
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where Q v and Q s represent source terms from volumetric and boundary contributions, 
respectively, k and C p denote the thermal conductivity and capacity of the solid material, 
respectively. In the case of simulating solid core with power generation, Q v depends on power 
generation distributions function employed. In the present study, three power distribution 
functions were examined: (i) pure sine function in the axial direction with pure cosine function 
in the radial direction; (ii) pure sine function in the axial direction with clipped cosine function 
in the radial direction; (iii) Clipped sine function in the axial direction with clipped cosine 
function in the radial direction. The temperature value at the fluid-solid interface is obtained by 
enforcing the heat flux continuity condition, i.e. Q s = -q w where q w is the heat flux from the fluid 
to the solid calculated by solving Eq. (3) including the turbulence effect if applicable. In order to 
achieve numerical stability and enforce the heat flux continuity condition, an implicit treatment 
of the temperature at the fluid-solid interface is employed. In this approach, Eq. (7) can be 
discretized as 



where the superscripts n+ 1 and n denote the values at the next and current time levels, 
respectively. 7* and T c are the temperatures at the fluid-solid interface and at the center of the 
solid cell next to the fluid-solid interface in respect. y„ represents the normal distance from the 
interface to the center of the solid cell next to the fluid-solid interface. The Vi factor on the left 


hand side of the above equation is because only half of the solid cell is involved in the control 
volume. It can be seen that this scheme can be applied to both transient and steady state 
simulations. For steady-state simulations, an acceleration factor can be used to improve 
convergence of heat conduction in the solid. Implementation of the implicit treatment has been 
validated 9 by comparing with the SINDA code. 


C, Porosity Model: 



In the present study, a porosity model was developed to represent the momentum and 
energy transport through an assembly of flow pipes, and the heat conduction through the solid 
material within a flow element. Hence, the porosity model will include separate temperatures 
and thermal conductivities for both the solid material and fluid flow. The momentum and energy 
equations for the fluid flow are the same as Eqs. 2 & 3, except the source terms will be modified 
to account for the extra friction loss and heat transfer. The source term of the momentum 
equation (Eq. 2) can be expressed as 
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where D is the drag force modeled by the porosity model, £ is the porosity factor for the porous 
region, V s is the superficial flow velocity through the porous region, d is diameter of the flow 
channel, and V is the total volume of the porous region. An empirical correlation of the friction 
coefficient (C/) for the flow through a pipe, known as the Blasius formula 24 , is used and 
expressed as follow. 

C f =0.079 IRe: 0 25 where Re. 

A# 

The source term for the energy equation (Eq. 3) can be calculated as 
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where Q is the heat sink/source due to the fluid-solid interaction in the porous region, Pr is the 
Prandtl number of the fluid, /is an empirical parameter will be tuned by comparing solutions of 
the porosity model and the detailed conjugate heat transfer model, and T s and T g are the 
temperatures of the solid and fluid at the same location, respectively. In addition, the source 
term of the heat conduction equation (Eq. 7) in the porous region can be calculated as 
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while Q s = 0. 


III. Numerical Meshes 

For the numerical simulations with the conjugate heat transfer model, two mesh systems 
need to be constructed. One of them is employed to model the entire computational domain 
including both the solid grain and the flow channel. The other one is used by the code to identify 
the region of solid material, and thus it can be very coarse. Hence, there is no need for a special 
grid generator to construct multi-block meshes for separating the solid and the fluid region, 
neither the identification for the fluid-solid interface. In the present study, various flow element 
geometries were simulated, such as different flow channel diameter, with and without coating 
layer in the solid core, and different flow channel length. Hence, a geometry/grid template, as 
shown in Figure 2, was developed based on MiniCAD 25 ' 27 to expedite the geometry and grid 
generation process by allowing the user to interactively change 1) the ratio of the gap between 
flow channels to the diameter of flow channel, 2) the ratio of the size of the flow element to the 
diameter of the flow channel, 3) thickness of the coating layer, and 4) the length of the flow 
element. The benefit of using MiniCAD is the ability to have graphical feedback while building 
the geometry and grid, and to easily move dimensional values from the geometry into the 



template to make it easy to create grids with modified values. In the present study, a 60° pi- 
section of a single flow element, which includes 3 and 1/6 flow channels, was simulated. 

MiniCAD is a graphical user interface (GUI) built on the Geometry Grid Toolkit 
(GGTK). The geometry that defines the dimensions of the solid grain was also exported in the 
stereolithography (STL) format, which can be used by an in-house grid generator to construct the 
hybrid mesh for the entire flow element. The process of building grid in MiniCAD involves 
three basic steps. Step one is to create or import the geometry. Step two is to pass the geometry 
to a grid topology to ensure that the mesh built will be geometrically watertight. Step 3 is to 
build the mesh based on the grid topology. Based on our investigation, we found the most 
efficient structured grid topology to model the solid grain is an O-type grid. In this case, the 
solid grain was decomposed into multiple blocks (5 blocks for cases without coating layer, and 
1 0 blocks for those with coating layer), where the O-grids of each block were used to model the 
solid grain and coating layer (if present) surrounding each flow channel as shown in Figure 2. A 
separate block is needed to identify the coating layer because it has different thermal properties 
from those of the solid grain. The boundary' of the O-grid geometry for each flow channel was 
built as a planar trimmed surface. This geometry was converted to the topology structure and 
then a 2-D grid was built on the faces. Since the grid for the solid grain is only used to identify 
the region of solid grain and is geometrically similar in the axial direction, it was only necessary 
to create the 3-D grid by extruding the 2-D grid in the axial direction, which can be seen in 
Figure 3. 



(a) 0. 1 mm coating layer (b) 0.2 mm coating layer 

Figure 2: Layout of the developed geometry/grid template (flow channel diameter = 2.05 mm) 



Figure 3: Structured grid system for the solid grain and coating layer of a flow element 

To accurately resolve the momentum and energy transport near the solid-fluid interface, a 
hybrid mesh topology was employed to generate the mesh for the entire flow element. The 
procedure used by our in-house grid generator is to construct the surface mesh before creating 
the 3-D volume mesh. For surface mesh generation, a direct advancing front method is used [28, 
29], A discrete model in STL format, imported from MiniCAD, is used as a background mesh, 
on which a new surface mesh suitable for computational simulations is created. Geometrical 
features are extracted based on a folding angle at each edge. A user specifies node distributions 
on the geometrical features, which form an initial front for the advancing front method. The 
surface triangulation is performed in the physical 3D space in order to check the quality of 
triangles easily. 

For volume mesh generation, two approaches were used. One is a method to create regular 
hybrid meshes based on an advancing layer method. In this method, a multiple marching 
direction approach is employed to create better control volumes around sharp convex comers 
[30]. Figure 4 shows a hybrid mesh inside a flow element, which has 335,207 nodes, 688,027 
tetrahedra, 21,476 prisms, 33,655 pyramids and 167,259 hexahedra. Quadrangles are used for 
the pipe surfaces. Multiple marching directions are defined at pipe inlets, which enable smooth 
transition of the mesh size there. 

The other volume mesh generation method is an extrusion method based on a 2-D mesh to create 
volume meshes more easily. Every cross-section across the symmetry axis is the same. First, a 
2-D mesh is created as shown in Figure 5. Second, it is extruded toward the symmetry axis, 
following the user-specified spacing. Third, boundary surfaces are numbered so that boundary 
conditions are set for CFD simulations. 




Figure 4: Regular hybrid meshes: (a) surface mesh; (b) cross-section of hybrid volume mesh 



Figure 5: 2D mesh for the extrusion method: solid part (yellow) and fluid part (green) 

IV. Results and Discussion 

The numerical results will be presented and discussed in the full length paper 
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